library(fpp2)
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
## ── Attaching packages ────────────────────────────────────────────── fpp2 2.4 ──
## ✓ ggplot2 3.3.3 ✓ fma 2.4
## ✓ forecast 8.14 ✓ expsmooth 2.3
##
library(xlsx)
library(tseries)
2. A classic example of a non-stationary series is the daily closing IBM stock price series (data set ibmclose). Use R to plot the daily closing prices for IBM stock and the ACF and PACF. Explain how each plot shows that the series is non-stationary and should be differenced.
ggtsdisplay(ibmclose)
# ACF plot shows that the autocorrelation values are bigger than critical value and decrease slowly. Also, r1 is large(near to 1) and positive.
# PACF plot shows that there is a strong correlation between IBM stock data and their 1 lagged values.
3. For the following series, find an appropriate Box-Cox transformation and order of differencing in order to obtain stationary data.
# a.
autoplot(usnetelec)
Box.test(diff(usnetelec), type = "Ljung-Box")
##
## Box-Ljung test
##
## data: diff(usnetelec)
## X-squared = 0.8508, df = 1, p-value = 0.3563
kpss.test(diff(usnetelec))
## Warning in kpss.test(diff(usnetelec)): p-value greater than printed p-value
##
## KPSS Test for Level Stationarity
##
## data: diff(usnetelec)
## KPSS Level = 0.15848, Truncation lag parameter = 3, p-value = 0.1
# kpss test result shows first differencing made the data stationary.
# b. usgdp
autoplot(usgdp)
Box.test(diff(usgdp), type = "Ljung-Box")
##
## Box-Ljung test
##
## data: diff(usgdp)
## X-squared = 39.187, df = 1, p-value = 3.85e-10
autoplot(diff(usgdp))
print("use ndiffs function to check the number of differencing.")
## [1] "use ndiffs function to check the number of differencing."
ndiffs(usgdp)
## [1] 2
autoplot(diff(diff(usgdp)))
Box.test(diff(diff(usgdp)), type = "Ljung-Box")
##
## Box-Ljung test
##
## data: diff(diff(usgdp))
## X-squared = 53.294, df = 1, p-value = 2.872e-13
ggAcf(diff(diff(usgdp)))
kpss.test(diff(diff(usnetelec)))
## Warning in kpss.test(diff(diff(usnetelec))): p-value greater than printed p-
## value
##
## KPSS Test for Level Stationarity
##
## data: diff(diff(usnetelec))
## KPSS Level = 0.098532, Truncation lag parameter = 3, p-value = 0.1
# c. mcopper
autoplot(mcopper)
lambda_mcopper <- BoxCox.lambda(mcopper)
autoplot(diff(BoxCox(mcopper, lambda_mcopper)))
Box.test(diff(BoxCox(mcopper, lambda_mcopper)),
type = "Ljung-Box")
##
## Box-Ljung test
##
## data: diff(BoxCox(mcopper, lambda_mcopper))
## X-squared = 57.517, df = 1, p-value = 3.353e-14
# Plot result looked like BoxCox transformation and first differencing made the data like white noise series.
ggAcf(diff(BoxCox(mcopper, lambda_mcopper)))
kpss.test(diff(BoxCox(mcopper, lambda_mcopper)))
## Warning in kpss.test(diff(BoxCox(mcopper, lambda_mcopper))): p-value greater
## than printed p-value
##
## KPSS Test for Level Stationarity
##
## data: diff(BoxCox(mcopper, lambda_mcopper))
## KPSS Level = 0.057275, Truncation lag parameter = 6, p-value = 0.1
# But kpss test result shows that differencing with Box-Cox transformation was enough to make the data stationary.
# Even if differencing with Box-Cox transformation didn't make the data like white noise series, it made the data stationary.
# d. enplanements
autoplot(enplanements)
# enplanements data have seasonality and increasing trend even if the number of enplanements fell in 2001. Therefore, I think that the data need seasonal differencing, too.
# The variations are bigger for bigger numbers.use Box-Cox transformation before differencing.
lambda_enplanements <- BoxCox.lambda(enplanements)
ndiffs(enplanements)
## [1] 1
nsdiffs(enplanements)
## [1] 1
# the data need 1 first differencing and 1 seasonal differencing.
autoplot(diff(diff(BoxCox(enplanements, lambda_enplanements),lag = 12)))
Box.test(diff(diff(BoxCox(enplanements, lambda_enplanements),lag = 12)),type = "Ljung-Box")
##
## Box-Ljung test
##
## data: diff(diff(BoxCox(enplanements, lambda_enplanements), lag = 12))
## X-squared = 29.562, df = 1, p-value = 5.417e-08
ggAcf(
diff(
diff(
BoxCox(enplanements, lambda_enplanements),
lag = 12
)
)
)
# There are still some autocorrelations left.
kpss.test(
diff(
diff(
BoxCox(enplanements, lambda_enplanements),
lag = 12
)
)
)
## Warning in kpss.test(diff(diff(BoxCox(enplanements, lambda_enplanements), : p-
## value greater than printed p-value
##
## KPSS Test for Level Stationarity
##
## data: diff(diff(BoxCox(enplanements, lambda_enplanements), lag = 12))
## KPSS Level = 0.042424, Truncation lag parameter = 5, p-value = 0.1
# But kpss test result shows that differencings with Box-Cox transformation was enough to make the data stationary. In enplanements data case, even if differencings with Box-Cox transformation didn't make the data like white noise series, it made the data stationary.
# e. visitors
autoplot(visitors)
# visitors data are similar to enplanements data. They have seasonality and increasing trend. It looked like they also need Box-Cox transformation, first and seasonal differencing.
lambda_visitors <- BoxCox.lambda(visitors)
ndiffs(visitors)
## [1] 1
nsdiffs(visitors)
## [1] 1
# visitors data need 1 first and 1 seasonal differencing.
autoplot(
diff(
diff(
BoxCox(visitors, lambda_visitors),
lag = 12
)
)
)
Box.test(
diff(
diff(
BoxCox(visitors, lambda_visitors),
lag = 12
)
),
type = "Ljung-Box"
)
##
## Box-Ljung test
##
## data: diff(diff(BoxCox(visitors, lambda_visitors), lag = 12))
## X-squared = 21.804, df = 1, p-value = 3.02e-06
# Plot result looked like BoxCox transformation and multiple differencings made the data like white noise series. But Ljung-Box test shows that it didn't.
ggAcf(
diff(
diff(
BoxCox(visitors, lambda_visitors),
lag = 12
)
)
)
# There are still some autocorrelations left.
kpss.test(
diff(
diff(
BoxCox(visitors, lambda_visitors),
lag = 12
)
)
)
## Warning in kpss.test(diff(diff(BoxCox(visitors, lambda_visitors), lag = 12))):
## p-value greater than printed p-value
##
## KPSS Test for Level Stationarity
##
## data: diff(diff(BoxCox(visitors, lambda_visitors), lag = 12))
## KPSS Level = 0.015833, Truncation lag parameter = 4, p-value = 0.1
# But kpss test result shows that differencings with Box-Cox transformation was enough to make the data stationary. In visitors data case, even if differencings with Box-Cox transformation didn't make the data like white noise series, it made the data stationary.
4. For the enplanements data, write down the differences you chose above using backshift operator notation.
the data needed 1 first difference, 1 seasonal difference after Box-Cox transformation. The model of the data can be written as ARIMA(0, 1, 0)(0, 1, 0)12 with Box-Cox transformation(lambda = -0.227).
The model expression using backshift operator notation B:
first equation : \(wt = (yt^(-0.227) - 1)/(-0.227)\)
second equation : \((1 - B)(1 - B^12)wt = et\), where et is a white noise series.
6. Use R to simulate and plot some data from simple ARIMA models.
# a. Use the following R code to generate data from an AR(1) model with phi1 = 0.6 and sigma^2 = 1. The process starts with y1 = 0.
y <- ts(numeric(100))
e <- rnorm(100)
for(i in 2:100){
y[i] <- 0.6*y[i-1] + e[i]
}
# b. Produce a time plot for the series. How does the plot change as you change phi1?
ar1generator <- function(phi1){
# generate 100 data points from an AR(1) model with input phi1.
y <- ts(numeric(100))
# error 'e's have variation sigma^2 as 1.
e <- rnorm(100)
for(i in 2:100){
y[i] <- phi1*y[i-1] + e[i]
}
return(y)
}
# produce plots changing phi1 value.
autoplot(ar1generator(0.3), series = "0.3") +
geom_line(size = 1, colour = "red") +
autolayer(y, series = "0.6", size = 1) +
autolayer(ar1generator(0.9), size = 1, series = "0.9") +
ylab("AR(1) models") +
guides(colour = guide_legend(title = "Phi1"))
# As phi increases, the variation of y increased.
# c. Write your own code to generate data from an MA(1) model with theta1 = 0.6 and sigma^2 = 1.
ma1generator <- function(theta1){
# generate 100 data points from an MA(1) model with input theta1.
y <- ts(numeric(100))
# error 'e's have variation sigma^2 as 1.
e <- rnorm(100)
for(i in 2:100){
y[i] <- theta1*e[i-1] + e[i]
}
return(y)
}
# d. Produce a time plot for the series. How does the plot change as you change theta1?
# produce plots changing theta1 value.
autoplot(ma1generator(0.3), series = "0.3") +
geom_line(size = 1, colour = "red") +
autolayer(y, series = "0.6", size = 1) +
autolayer(ar1generator(0.9), size = 1, series = "0.9") +
ylab("MA(1) models") +
guides(colour = guide_legend(title = "Theta1"))
# As theta increases, the variation of y increased.
# e. Generate data from an ARMA(1,1) model with phi1 = 0.6, theta1 = 0.6 and sigma^2 = 1.
y_arima.1.0.1 <- ts(numeric(50))
e <- rnorm(50)
for(i in 2:50){
y_arima.1.0.1[i] <- 0.6*y_arima.1.0.1[i-1] + 0.6*e[i-1] + e[i]
}
# f. Generate data from an AR(2) model with phi1 = -0.8, phi2 = 0.3 and sigma^2 = 1. (Note that these parameters will give a non-stationary series.)
y_arima.2.0.0 <- ts(numeric(50))
e <- rnorm(50)
for(i in 3:50){
y_arima.2.0.0[i] <- -0.8*y_arima.2.0.0[i-1] + 0.3*y_arima.2.0.0[i-2] + e[i]
}
# g. Graph the latter two series and compare them.
autoplot(y_arima.1.0.1, series = "ARMA(1, 1)") +
autolayer(y_arima.2.0.0, series = "AR(2)") +
ylab("y") +
guides(colour = guide_legend(title = "Models"))
autoplot(y_arima.1.0.1)
# data from an AR(2) model increased with oscillation. They are non-staionary data. But data from an ARMA(1, 1) model were stationary.
7. Consider the number of women murdered each year (per 100,000 standard population) in the United States. (Data set wmurders).
# a. By studying appropriate graphs of the series in R, find an appropriate ARIMA(p,d,q) model for these data.
autoplot(wmurders)
# It looked like the data don't need seasonal differencing or Box-Cox transformation.
autoplot(diff(wmurders))
# It looked like 1 more differencing would be needed to make the data stationary. Differenced data slowly go to minus infinity.
ndiffs(wmurders)
## [1] 2
# ndiffs function shows that the data need 2 differencing.
autoplot(diff(wmurders, differences = 2))
kpss.test(diff(wmurders, differences = 2))
## Warning in kpss.test(diff(wmurders, differences = 2)): p-value greater than
## printed p-value
##
## KPSS Test for Level Stationarity
##
## data: diff(wmurders, differences = 2)
## KPSS Level = 0.045793, Truncation lag parameter = 3, p-value = 0.1
# twice differencing made the data stationary.
diff(wmurders, differences = 2) %>% ggtsdisplay()
# b. Should you include a constant in the model? Explain.
# ARIMA model of the data includes twice differencing. If there is a constant in the model, twice integrated contant will yield quadratic trend, which is dangerous for forecasting. Therefore I won't include a constant in the model.
# c. Write this model in terms of the backshift operator.
# (1 - B)^2*yt = (1 + theta1*B + theta2*B^2)*et
# d. Fit the model using R and examine the residuals. Is the model satisfactory?
wmurders_arima.0.2.2 <- Arima(wmurders,
order = c(0, 2, 2))
checkresiduals(wmurders_arima.0.2.2)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,2,2)
## Q* = 11.764, df = 8, p-value = 0.1621
##
## Model df: 2. Total lags used: 10
# The residuals of the model can be thought of as white noise series. A little sorry that they aren't normally distributed. But it is satisfactory to get them.
# e. Forecast three times ahead. Check your forecasts by hand to make sure that you know how they have been calculated.
fc_wmurders_arima.0.2.2 <- forecast(
wmurders_arima.0.2.2, h = 3
)
# forecasts by Arima function
fc_wmurders_arima.0.2.2$mean
## Time Series:
## Start = 2005
## End = 2007
## Frequency = 1
## [1] 2.480525 2.374890 2.269256
# get forecasts by manual calculation
fc_wmurders_arima.0.2.2$model
## Series: wmurders
## ARIMA(0,2,2)
##
## Coefficients:
## ma1 ma2
## -1.0181 0.1470
## s.e. 0.1220 0.1156
##
## sigma^2 estimated as 0.04702: log likelihood=6.03
## AIC=-6.06 AICc=-5.57 BIC=-0.15
years <- length(wmurders)
e <- fc_wmurders_arima.0.2.2$residuals
fc1 <- 2*wmurders[years] - wmurders[years - 1] - 1.0181*e[years] + 0.1470*e[years - 1]
fc2 <- 2*fc1 - wmurders[years] + 0.1470*e[years]
fc3 <- 2*fc2 - fc1
# forecasts by manual calculation
c(fc1, fc2, fc3)
## [1] 2.480523 2.374887 2.269252
# the forecasts are almost similar to the ones got by Arima function.
# f. Create a plot of the series with forecasts and prediction intervals for the next three periods shown.
autoplot(fc_wmurders_arima.0.2.2)
# g. Does auto.arima give the same model you have chosen? If not, which model do you think is better?
fc_wmurders_autoarima <- forecast(
auto.arima(wmurders), h = 3
)
# Without RMSE, all errors show that ARIMA(0, 2, 2) is better than ARIMA(1, 2, 1).
accuracy(fc_wmurders_arima.0.2.2)
## ME RMSE MAE MPE MAPE MASE
## Training set -0.0113461 0.2088162 0.1525773 -0.2403396 4.331729 0.9382785
## ACF1
## Training set -0.05094066
accuracy(fc_wmurders_autoarima)
## ME RMSE MAE MPE MAPE MASE
## Training set -0.01065956 0.2072523 0.1528734 -0.2149476 4.335214 0.9400996
## ACF1
## Training set 0.02176343
# try using auto.arima function with stepwise and approximation options false.
fc_wmurders_autoarima2 <- forecast(
auto.arima(wmurders, stepwise = FALSE, approximation = FALSE),
h = 3
)
# It is ARIMA(0, 2, 3) model.
accuracy(fc_wmurders_autoarima2)
## ME RMSE MAE MPE MAPE MASE
## Training set -0.01336585 0.2016929 0.1531053 -0.3332051 4.387024 0.9415259
## ACF1
## Training set -0.03193856
# In this case, some errors were better while others were worse. I'll check residuals and ACF, PACF plots.
ggtsdisplay(diff(wmurders, differences = 2))
# It looked like that the data are similar to ARIMA(0, 2, 2) rather than ARIMA(0, 2, 3).
checkresiduals(fc_wmurders_arima.0.2.2)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,2,2)
## Q* = 11.764, df = 8, p-value = 0.1621
##
## Model df: 2. Total lags used: 10
checkresiduals(fc_wmurders_autoarima2)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,2,3)
## Q* = 10.706, df = 7, p-value = 0.152
##
## Model df: 3. Total lags used: 10
# almost similar residuals.
# Therefore I'll choose ARIMA(0, 2, 2).
8. Consider the total international visitors to Australia (in millions) for the period 1980-2015. (Data set austa.)
# a. Use auto.arima to find an appropriate ARIMA model. What model was selected. Check that the residuals look like white noise. Plot forecasts for the next 10 periods.
autoplot(austa)
fc_austa_autoarima <- forecast(
auto.arima(austa), h = 10
)
fc_austa_autoarima$model
## Series: austa
## ARIMA(0,1,1) with drift
##
## Coefficients:
## ma1 drift
## 0.3006 0.1735
## s.e. 0.1647 0.0390
##
## sigma^2 estimated as 0.03376: log likelihood=10.62
## AIC=-15.24 AICc=-14.46 BIC=-10.57
# ARIMA(0, 1, 1) with drift model was chosen.
checkresiduals(fc_austa_autoarima)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,1,1) with drift
## Q* = 2.297, df = 5, p-value = 0.8067
##
## Model df: 2. Total lags used: 7
# The residuals are like white noise.
autoplot(fc_austa_autoarima)
# b. Plot forecasts from an ARIMA(0,1,1) model with no drift and compare these to part (a). Remove the MA term and plot again.
fc_austa_arima.0.1.1 <- forecast(
Arima(austa, order = c(0, 1, 1)), h = 10
)
autoplot(fc_austa_arima.0.1.1)
fc_austa_arima.0.1.0 <- forecast(
Arima(austa, order = c(0, 1, 0)), h = 10
)
autoplot(fc_austa_arima.0.1.0)
# the forecasts of both models are like the result of naive forecast. Increasing trend isn't reflected in the forecasts.
fc_austa_arima.0.1.1$upper - fc_austa_arima.0.1.0$upper
## Time Series:
## Start = 2016
## End = 2025
## Frequency = 1
## fc_austa_arima.0.1.1$upper.80% fc_austa_arima.0.1.1$upper.95%
## 2016 0.1138117 0.09101057
## 2017 0.2039418 0.22885273
## 2018 0.2527212 0.30345426
## 2019 0.2886706 0.35843413
## 2020 0.3180942 0.40343373
## 2021 0.3434784 0.44225553
## 2022 0.3660754 0.47681466
## 2023 0.3866116 0.50822207
## 2024 0.4055495 0.53718500
## 2025 0.4232034 0.56418442
fc_austa_arima.0.1.0$lower - fc_austa_arima.0.1.1$lower
## Time Series:
## Start = 2016
## End = 2025
## Frequency = 1
## fc_austa_arima.0.1.0$lower.80% fc_austa_arima.0.1.0$lower.95%
## 2016 -0.199956384 -0.22275751
## 2017 -0.109826244 -0.08491535
## 2018 -0.061046922 -0.01031382
## 2019 -0.025097519 0.04466605
## 2020 0.004326137 0.08966565
## 2021 0.029710347 0.12848745
## 2022 0.052307350 0.16304658
## 2023 0.072843552 0.19445399
## 2024 0.091781392 0.22341692
## 2025 0.109435361 0.25041634
# But prediction interval of ARIMA(0, 1, 1) model was generally larger than the one of ARIMA(0, 1, 0) model. I think that it is because of one more error term in ARIMA(0, 1, 1) model.
# c. Plot forecasts from an ARIMA(2,1,3) model with drift. Remove the constant and see what happens.
fc_austa_arima.2.1.3.drift <- forecast(
Arima(austa, order = c(2, 1, 3), include.drift = TRUE),
h = 10
)
autoplot(fc_austa_arima.2.1.3.drift)
# The forecasts are increasing, but the speed of the increase is decreasing.
drift_austa <- fc_austa_arima.2.1.3.drift$model$coef[6]
fc_austa_arima.2.1.3.nodrift <- fc_austa_arima.2.1.3.drift$mean - drift_austa*seq_len(10)
autoplot(fc_austa_arima.2.1.3.drift) +
autolayer(fc_austa_arima.2.1.3.nodrift)
# Without drift constant, the forecasts are unlikely.
# d. Plot forecasts from an ARIMA(0,0,1) model with a constant. Remove the MA term and plot again.
fc_austa_arima.0.0.1.const <- forecast(
Arima(
austa, order = c(0, 0, 1), include.constant = TRUE
),
h = 10
)
autoplot(fc_austa_arima.0.0.1.const)
# the forecasts are fastly decreased to the mean of the data history.
fc_austa_arima.0.0.0.const <- forecast(
Arima(austa, order = c(0, 0, 0), include.constant = TRUE),
h = 10
)
autoplot(fc_austa_arima.0.0.0.const)
# All of the forecasts are the mean of the data history. It is like the result of mean method.
# e. Plot forecasts from an ARIMA(0,2,1) model with no constant.
fc_austa_arima.0.2.1 <- forecast(
Arima(austa, order = c(0, 2, 1)),
h = 10
)
autoplot(fc_austa_arima.0.2.1)
# the forecasts show increasing trend. PI is being larger for the farther future forecast.
9. For the usgdp series:
# a. if necessary, find a suitable Box-Cox transformation for the data;
autoplot(usgdp)
autoplot(BoxCox(usgdp, BoxCox.lambda(usgdp)))
lambda_usgdp <- BoxCox.lambda(usgdp)
# b.fit a suitable ARIMA model to the transformed data using auto.arima();
usgdp_autoarima <- auto.arima(usgdp,
lambda = lambda_usgdp)
autoplot(usgdp, series = "Data") +
autolayer(usgdp_autoarima$fitted, series = "Fitted")
# It looked like the model fits well to the data.
usgdp_autoarima
## Series: usgdp
## ARIMA(2,1,0) with drift
## Box Cox transformation: lambda= 0.366352
##
## Coefficients:
## ar1 ar2 drift
## 0.2795 0.1208 0.1829
## s.e. 0.0647 0.0648 0.0202
##
## sigma^2 estimated as 0.03518: log likelihood=61.56
## AIC=-115.11 AICc=-114.94 BIC=-101.26
#ARIMA(2, 1, 0) with drift model after Box-Cox transformation.
# c. try some other plausible models by experimenting with the orders chosen;
ndiffs(BoxCox(usgdp, lambda_usgdp))
## [1] 1
# the data need 1 first differencing to be stationary.
ggtsdisplay(diff(BoxCox(usgdp, lambda_usgdp)))
# ACF plot shows sinusoidal decrease while PACF plot shows significant spikes at lag 1 and 12. I think that I can ignore the spike at lag 12 because the data are aggregated quarterly, not monthly. Therefore, I'll experiment with ARIMA(1, 1, 0) model.
usgdp_arima.1.1.0 <- Arima(
usgdp, lambda = lambda_usgdp, order = c(1, 1, 0)
)
usgdp_arima.1.1.0
## Series: usgdp
## ARIMA(1,1,0)
## Box Cox transformation: lambda= 0.366352
##
## Coefficients:
## ar1
## 0.6326
## s.e. 0.0504
##
## sigma^2 estimated as 0.04384: log likelihood=34.39
## AIC=-64.78 AICc=-64.73 BIC=-57.85
autoplot(usgdp, series = "Data") +
autolayer(usgdp_arima.1.1.0$fitted, series = "Fitted")
# I'll also try ARIMA(1, 1, 0) with drift model.
usgdp_arima.1.1.0.drift <- Arima(
usgdp, lambda = lambda_usgdp, order = c(1, 1, 0),
include.drift = TRUE
)
usgdp_arima.1.1.0.drift
## Series: usgdp
## ARIMA(1,1,0) with drift
## Box Cox transformation: lambda= 0.366352
##
## Coefficients:
## ar1 drift
## 0.3180 0.1831
## s.e. 0.0619 0.0179
##
## sigma^2 estimated as 0.03555: log likelihood=59.83
## AIC=-113.66 AICc=-113.56 BIC=-103.27
autoplot(usgdp, series = "Data") +
autolayer(usgdp_arima.1.1.0.drift$fitted, series = "Fitted")
# It looked like that these models also fit well to the data.
# d. choose what you think is the best model and check the residual diagnostics;
accuracy(usgdp_autoarima)
## ME RMSE MAE MPE MAPE MASE
## Training set 1.195275 39.2224 29.29521 -0.01363259 0.6863491 0.1655687
## ACF1
## Training set -0.03824844
accuracy(usgdp_arima.1.1.0)
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set 15.45449 45.49569 35.08393 0.3101283 0.7815664 0.198285 -0.3381619
accuracy(usgdp_arima.1.1.0.drift)
## ME RMSE MAE MPE MAPE MASE
## Training set 1.315796 39.90012 29.5802 -0.01678591 0.6834509 0.1671794
## ACF1
## Training set -0.08544569
# Some errors show that ARIMA(2, 1, 0) with drift is the best model while others show that ARIMA(1, 1, 0) with drift is the best. Check the residuals of both cases.
checkresiduals(usgdp_autoarima)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(2,1,0) with drift
## Q* = 6.5772, df = 5, p-value = 0.254
##
## Model df: 3. Total lags used: 8
checkresiduals(usgdp_arima.1.1.0.drift)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,1,0) with drift
## Q* = 10.274, df = 6, p-value = 0.1136
##
## Model df: 2. Total lags used: 8
# In either case, the residuals are like white noise series and are not normally distributed.
# With the model, RMSE and MASE values were lower. And there wasn't significant spike at lag 2 in ACF plot of ARIMA(2, 1, 0) with drift model, even if it exists in ARIMA(1, 1, 0) with drift model.
# e. produce forecasts of your fitted model. Do the forecasts look reasonable?
fc_usgdp_autoarima <- forecast(
usgdp_autoarima
)
autoplot(fc_usgdp_autoarima)
# It looked like the forecasts are reasonable.
# f. compare the results with what you would obtain using ets() (with no transformation).
fc_usgdp_ets <- forecast(
ets(usgdp)
)
autoplot(fc_usgdp_ets)
# It looked like these forecasts are more likely than the ones with ARIMA model. When trend is obvious, is ETS better than ARIMA model? I wonder about it.
10. Consider austourists, the quarterly number of international tourists to Australia for the period 1999-2010. (Data set austourists.)
# a. Describe the time plot.
autoplot(austourists)
# the data have strong seasonality and increasing trend. Also the size of variations increased as the number increased.
# b. What can you learn from the ACF graph?
ggAcf(austourists)
# autocorrelations are slowly decreasing. And the values at the lags of multiple of 4 were big compared to the others.
# c. What can you learn from the PACF graph?
ggPacf(austourists)
# there are 5 significant spikes, and then no significant spikes thereafter (apart from one at lag 8, which are probably related with quarterly seasonality).
# d. Produce plots of the seasonally differenced data (1 - B^4)Yt. What model do these graphs suggest?
ggtsdisplay(diff(austourists, lag = 4))
# the seasonally differenced data are looked like to need at least one more differencing to make it stationary.
# For the values at the lags of multiple of 4, there are just significant spikes at lag 4. It is same for ACF and PACF plots. The order of seasonal ARIMA model can be (1, 1, 0)[4] or (0, 1, 1)[4]. I'll choose (1, 1, 0)[4] order.
# Disregarding the values at lag 4, autocorrelation values are looked like decreasing sinusoidally while partial autocorrelation values have spikes at lag 1 and 5. There aren't significant spikes at lag 2 and 3.
# I think that the spike at lag 5 doesn't mean that there are still important information unused at 5th lagged values. It should've been positively significant because of the big negatively significant spike at lag 4 after seasonal differencing.
# I can find this by drawing PACF plot after doing one more differencing.
ggtsdisplay(diff(diff(austourists, lag = 4)))
# Therefore I suggest ARIMA(1, 1, 0)(1, 1, 0)[4] model.
# e. Does auto.arima give the same model that you chose? If not, which model do you think is better?
fc_austourists_autoarima <- forecast(
auto.arima(austourists)
)
fc_austourists_autoarima$model
## Series: austourists
## ARIMA(1,0,0)(1,1,0)[4] with drift
##
## Coefficients:
## ar1 sar1 drift
## 0.4705 -0.5305 0.5489
## s.e. 0.1154 0.1122 0.0864
##
## sigma^2 estimated as 5.15: log likelihood=-142.48
## AIC=292.97 AICc=293.65 BIC=301.6
# auto.arima gave ARIMA(1, 0, 0)(1, 1, 0)[4] model.
fc_austourists_arima.1.1.0.1.1.0.4 <- forecast(
Arima(austourists,
order = c(1, 1, 0),
seasonal = c(1, 1, 0))
)
autoplot(fc_austourists_autoarima)
autoplot(fc_austourists_arima.1.1.0.1.1.0.4)
# ARIMA(1, 1, 0)(1, 1, 0)[4] shows more fastly increasing trend.
accuracy(fc_austourists_autoarima)
## ME RMSE MAE MPE MAPE MASE
## Training set 0.02200144 2.149384 1.620917 -0.7072593 4.388288 0.5378929
## ACF1
## Training set -0.06393238
accuracy(fc_austourists_arima.1.1.0.1.1.0.4)
## ME RMSE MAE MPE MAPE MASE
## Training set 0.07547467 2.288429 1.661709 -0.3718168 4.49153 0.5514295
## ACF1
## Training set -0.03831574
# ARIMA(1, 0, 0)(1, 1, 0)[4] with drift model was fitted better. Therefore I think that this model is better than ARIMA(1, 1, 0)(1, 1, 0)[4] model.
checkresiduals(fc_austourists_autoarima)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,0,0)(1,1,0)[4] with drift
## Q* = 4.0937, df = 5, p-value = 0.536
##
## Model df: 3. Total lags used: 8
# The residuals are like white noise series.
# f. Write the model in terms of the backshift operator, then without using the backshift operator.
fc_austourists_autoarima$model
## Series: austourists
## ARIMA(1,0,0)(1,1,0)[4] with drift
##
## Coefficients:
## ar1 sar1 drift
## 0.4705 -0.5305 0.5489
## s.e. 0.1154 0.1122 0.0864
##
## sigma^2 estimated as 5.15: log likelihood=-142.48
## AIC=292.97 AICc=293.65 BIC=301.6
# (1 - phi1*B)(1 - phis1*B)(1 - B^4)(yt - c*t) = et
# c = drift*(1 - phi1)(1 - phis1)*m^D = 1.7793
# (1 - phi1*B - phis1*B + phi1*phis1*B^2)(1 - B^4)(yt - c*t) =
# (1 - phi1*B - phis1*B + phi1*phis1*B^2 - B^4 + phi1*B^5 + phis1*B^5 - phi1*phis1*B^6)(yt - c*t) = et
# yt = c + (phi1 + phis1)*yt-1 - phi1*phis1*yt-2 + yt-4 - (phi1 + phis1)*yt-5 + phi1*phis1*yt-6 + et
# yt = 1.7793 - 0.06*yt-1 + 0.2496*yt-2 + yt-4 + 0.06*yt-5 - 0.2496*yt-6 + et
11. Consider the total net generation of electricity (in billion kilowatt hours) by the U.S. electric industry (monthly for the period January 1973 - June 2013). (Data set usmelec.) In general there are two peaks per year: in mid-summer and mid-winter.
# a. Examine the 12-month moving average of this series to see what kind of trend is involved.
usmelec_ma2x12 <- ma(usmelec, order = 12, centre = TRUE)
autoplot(usmelec, series = "Data") +
autolayer(usmelec_ma2x12, series = "2X12-MA") +
ylab(expression(paste("Electricity(x", 10^{9}, "KWh)"))) +
ggtitle("Monthly total net generation of electricity") +
scale_color_discrete(breaks = c("Data", "2X12-MA"))
## Warning: Removed 12 row(s) containing missing values (geom_path).
# Total net generation amount increased first but stoped increasing from about 2008.
# b. Do the data need transforming? If so, find a suitable transformation.
# The data show bigger variation for bigger amount. Therefore I think that Box-Cox transformation would be suitable for the data.
lambda_usmelec <- BoxCox.lambda(usmelec)
# c. Are the data stationary? If not, find an appropriate differencing which yields stationary data.
# The data are non-stationary.
ndiffs(usmelec)
## [1] 1
nsdiffs(usmelec)
## [1] 1
# I need to do 1 seasonal differencing to make the data stationary. If seasonal differencing isn't enough to make them stationary, I need to do first differencing, too.
# d. Identify a couple of ARIMA models that might be useful in describing the time series. Which of your models is the best according to their AIC values?
ggtsdisplay(diff(
BoxCox(usmelec, lambda_usmelec),
lag = 12
))
# Definitely, I need to use first differencing, too.
ggtsdisplay(
diff(
diff(
BoxCox(usmelec, lambda_usmelec),
lag = 12
)
)
)
# I think that ARIMA(0, 1, 2)(0, 1, 1)[12] with Box-Cox transformation model might describe the data well. I'll try ARIMA(0, 1, 3)(0, 1, 1)[12] with Box-Cox transformation model, too.
usmelec_arima.0.1.2.0.1.1.12 <- Arima(
usmelec,
lambda = lambda_usmelec,
order = c(0, 1, 2),
seasonal = c(0, 1, 1)
)
usmelec_arima.0.1.3.0.1.1.12 <- Arima(
usmelec,
lambda = lambda_usmelec,
order = c(0, 1, 3),
seasonal = c(0, 1, 1)
)
usmelec_arima.0.1.2.0.1.1.12$aic
## [1] -5081.506
usmelec_arima.0.1.3.0.1.1.12$aic
## [1] -5080.441
# ARIMA(0, 1, 2)(0, 1, 1)[12] with Box-Cox transformation model was the best.
# e. Estimate the parameters of your best model and do diagnostic testing on the residuals. Do the residuals resemble white noise? If not, try to find another ARIMA model which fits better.
usmelec_arima.0.1.2.0.1.1.12
## Series: usmelec
## ARIMA(0,1,2)(0,1,1)[12]
## Box Cox transformation: lambda= -0.5738331
##
## Coefficients:
## ma1 ma2 sma1
## -0.4317 -0.2552 -0.8536
## s.e. 0.0439 0.0440 0.0261
##
## sigma^2 estimated as 1.284e-06: log likelihood=2544.75
## AIC=-5081.51 AICc=-5081.42 BIC=-5064.87
#theta1 = -0.4317, theta2 = -0.2552, phis1 = -0.8536
checkresiduals(usmelec_arima.0.1.2.0.1.1.12)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,1,2)(0,1,1)[12]
## Q* = 32.525, df = 21, p-value = 0.05176
##
## Model df: 3. Total lags used: 24
# Ljung-Box test result shows that the residuals can be thought of as white noise. And they are normally distributed.
# I want to know what model was selected if I used auto.arima function. I'll try it.
usmelec_autoarima <- auto.arima(
usmelec,
lambda = lambda_usmelec
)
usmelec_autoarima
## Series: usmelec
## ARIMA(1,1,3)(2,1,1)[12]
## Box Cox transformation: lambda= -0.5738331
##
## Coefficients:
## ar1 ma1 ma2 ma3 sar1 sar2 sma1
## 0.3952 -0.8194 -0.0468 0.0414 0.0403 -0.0934 -0.8462
## s.e. 0.3716 0.3734 0.1707 0.1079 0.0560 0.0533 0.0343
##
## sigma^2 estimated as 1.278e-06: log likelihood=2547.75
## AIC=-5079.5 AICc=-5079.19 BIC=-5046.23
# The result is ARIMA(2, 1, 4)(0, 0, 2)[12] with drift after Box-Cox transformation model. AIC is -4722. But I can't compare the AIC value with what I got above, because the number of differencing was different(Differencing changes the way the likelihood is computed).
checkresiduals(usmelec_autoarima)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,1,3)(2,1,1)[12]
## Q* = 25.995, df = 17, p-value = 0.07454
##
## Model df: 7. Total lags used: 24
# And the residuals aren't like white noise. Therefore I'll choose ARIMA(0, 1, 2)(0, 1, 1)[12] with Box-Cox transformation model.
# g. How many years of forecasts do you think are sufficiently accurate to be usable?
# In usmelec data case, even 4 years of forecasts were sufficiently accurate to be usable. I think that it happened because the pattern in the data almost didn't change.
12. For the mcopper data:
# a. if necessary, find a suitable Box-Cox transformation for the data;
autoplot(mcopper)
# they are monthly data but there isn't seasonality in them.
autoplot(BoxCox(mcopper, BoxCox.lambda(mcopper)))
# It looked like Box-Cox transformation makes the variations in the data evenly over time. Therefore I'm going to use the transformation.
lambda_mcopper <- BoxCox.lambda(mcopper)
# b. fit a suitable ARIMA model to the transformed data using auto.arima();
mcopper_autoarima <- auto.arima(
mcopper,
lambda = lambda_mcopper
)
mcopper_autoarima
## Series: mcopper
## ARIMA(0,1,1)
## Box Cox transformation: lambda= 0.1919047
##
## Coefficients:
## ma1
## 0.3720
## s.e. 0.0388
##
## sigma^2 estimated as 0.04997: log likelihood=45.05
## AIC=-86.1 AICc=-86.08 BIC=-77.43
# auto.arima yielded ARIMA(0, 1, 1) with Box-Cox transformation model. AICc was -86.08.
# c. try some other plausible models by experimenting with the orders chosen;
ndiffs(mcopper)
## [1] 1
nsdiffs(mcopper)
## [1] 0
# the data need 1 first differencing.
ggtsdisplay(diff(mcopper))
# It looked like autocorrelation values are sinusoidally decreasing. I'll choose ARIMA model's order as (1, 1, 0) and (5, 1, 0).
mcopper_arima.1.1.0 <- Arima(
mcopper, order = c(1, 1, 0), lambda = lambda_mcopper
)
mcopper_arima.1.1.0
## Series: mcopper
## ARIMA(1,1,0)
## Box Cox transformation: lambda= 0.1919047
##
## Coefficients:
## ar1
## 0.3231
## s.e. 0.0399
##
## sigma^2 estimated as 0.05091: log likelihood=39.83
## AIC=-75.66 AICc=-75.64 BIC=-66.99
# AICc was -75.64.
mcopper_arima.5.1.0 <- Arima(
mcopper, order = c(5, 1, 0), lambda = lambda_mcopper
)
mcopper_arima.5.1.0
## Series: mcopper
## ARIMA(5,1,0)
## Box Cox transformation: lambda= 0.1919047
##
## Coefficients:
## ar1 ar2 ar3 ar4 ar5
## 0.3706 -0.1475 0.0609 -0.0038 0.0134
## s.e. 0.0421 0.0450 0.0454 0.0450 0.0422
##
## sigma^2 estimated as 0.05028: log likelihood=45.32
## AIC=-78.63 AICc=-78.48 BIC=-52.63
# AICc was -78.48.
# I'll try auto.arima function without approximation and stepwise options.
mcopper_autoarima2 <- auto.arima(
mcopper, lambda = lambda_mcopper,
approximation = FALSE, stepwise = FALSE
)
mcopper_autoarima2
## Series: mcopper
## ARIMA(0,1,1)
## Box Cox transformation: lambda= 0.1919047
##
## Coefficients:
## ma1
## 0.3720
## s.e. 0.0388
##
## sigma^2 estimated as 0.04997: log likelihood=45.05
## AIC=-86.1 AICc=-86.08 BIC=-77.43
# the result model is the same as when I didn't use the options.
# d. choose what you think is the best model and check the residual diagnostics;
# When I compared AICc values, I got the smallest when I used auto.arima function. (I could've used AICc in comparing because the differencing was the same for all models I chose.) Best model is ARIMA(0, 1, 1) with Box-Cox transformation.
checkresiduals(mcopper_autoarima)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,1,1)
## Q* = 22.913, df = 23, p-value = 0.4659
##
## Model df: 1. Total lags used: 24
# The residuals are like white noise. I'll select the model.
# e. produce forecasts of your fitted model. Do the forecasts look reasonable?
fc_mcopper_autoarima <- forecast(
mcopper_autoarima
)
autoplot(fc_mcopper_autoarima)
# The forecasts aren't reasonable.
# I'll try other models I made.
fc_mcopper_arima.1.1.0 <- forecast(
mcopper_arima.1.1.0
)
autoplot(fc_mcopper_arima.1.1.0)
# got almost same result
fc_mcopper_arima.5.1.0 <- forecast(
mcopper_arima.5.1.0
)
autoplot(fc_mcopper_arima.5.1.0)
# got almost same result, too
# f. compare the results with what you would obtain using ets() (with no transformation).
fc_mcopper_ets <- forecast(
ets(mcopper)
)
autoplot(fc_mcopper_ets)
# These forecasts are more reasonable than what I got above.
13. Choose one of the following seasonal time series: hsales, auscafe, qauselec, qcement, qgas.
# a. Do the data need transforming? If so, find a suitable transformation.
# I'll analyze and forecast qauselec data.
autoplot(qauselec)
# The data need Box-Cox transformation to make the variations evenly over time.
lambda_qauselec <- BoxCox.lambda(qauselec)
# b. Are the data stationary? If not, find an appropriate differencing which yields stationary data.
# The data have strong seasonality and increasing trend. It means that the data aren't stationary.
nsdiffs(qauselec)
## [1] 1
ndiffs(qauselec)
## [1] 1
# The data need 1 seasonal differencing.
kpss.test(diff(qauselec, lag = 4))
##
## KPSS Test for Level Stationarity
##
## data: diff(qauselec, lag = 4)
## KPSS Level = 0.39382, Truncation lag parameter = 4, p-value = 0.07982
# The data don't need first differencing. But I'll try with first differencing, too.
# c. Identify a couple of ARIMA models that might be useful in describing the time series. Which of your models is the best according to their AIC values?
ggtsdisplay(diff(
BoxCox(qauselec, lambda_qauselec), lag = 4
))
ggtsdisplay(diff(diff(
BoxCox(qauselec, lambda_qauselec), lag = 4
)))
# suggest models:
# ARIMA(0, 0, 1)(0, 1, 1)[4]
qauselec_arima0.0.1.0.1.1.4 <- Arima(
qauselec, lambda = lambda_qauselec,
order = c(0, 0, 1), seasonal = c(0, 1, 1)
)
qauselec_arima0.0.1.0.1.1.4
## Series: qauselec
## ARIMA(0,0,1)(0,1,1)[4]
## Box Cox transformation: lambda= 0.5195274
##
## Coefficients:
## ma1 sma1
## 0.6495 0.2738
## s.e. 0.0746 0.0642
##
## sigma^2 estimated as 0.03639: log likelihood=51.5
## AIC=-97.01 AICc=-96.89 BIC=-86.91
# AIC = -97.01
# ARIMA(0, 1, 1)(0, 1, 1)[4]
qauselec_arima0.1.1.0.1.1.4 <- Arima(
qauselec, lambda = lambda_qauselec,
order = c(0, 1, 1), seasonal = c(0, 1, 1)
)
qauselec_arima0.1.1.0.1.1.4
## Series: qauselec
## ARIMA(0,1,1)(0,1,1)[4]
## Box Cox transformation: lambda= 0.5195274
##
## Coefficients:
## ma1 sma1
## -0.4976 -0.6910
## s.e. 0.0772 0.0418
##
## sigma^2 estimated as 0.01435: log likelihood=149.29
## AIC=-292.59 AICc=-292.47 BIC=-282.5
#AIC = -292.59
# ARIMA(0, 1, 1)(0, 1, 2)[4]
qauselec_arima0.1.1.0.1.2.4 <- Arima(
qauselec, lambda = lambda_qauselec,
order = c(0, 1, 1), seasonal = c(0, 1, 2)
)
qauselec_arima0.1.1.0.1.2.4
## Series: qauselec
## ARIMA(0,1,1)(0,1,2)[4]
## Box Cox transformation: lambda= 0.5195274
##
## Coefficients:
## ma1 sma1 sma2
## -0.5037 -0.7993 0.1249
## s.e. 0.0710 0.0870 0.0864
##
## sigma^2 estimated as 0.01425: log likelihood=150.36
## AIC=-292.73 AICc=-292.53 BIC=-279.28
#AIC = -292.73
# try using auto.arima function
qauselec_autoarima <- auto.arima(
qauselec, lambda = lambda_qauselec
)
qauselec_autoarima
## Series: qauselec
## ARIMA(1,1,1)(1,1,2)[4]
## Box Cox transformation: lambda= 0.5195274
##
## Coefficients:
## ar1 ma1 sar1 sma1 sma2
## 0.2523 -0.6905 0.8878 -1.6954 0.7641
## s.e. 0.1562 0.1279 0.0864 0.0983 0.0772
##
## sigma^2 estimated as 0.01345: log likelihood=156.42
## AIC=-300.84 AICc=-300.43 BIC=-280.67
# AIC = -300.84
# According to AIC values, ARIMA(1, 1, 1)(1, 1, 2)[4] with Box-Cox transformation is the best model. But above 2 models and the model made by auto.arima function used different number of differencings.
# d. Estimate the parameters of your best model and do diagnostic testing on the residuals. Do the residuals resemble white noise? If not, try to find another ARIMA model which fits better.
# phi1 = 0.2523, theta1 = -0.6905, phis1 = 0.8878, thetas1 = -1.6954, thetas2 = -0.7641
checkresiduals(qauselec_autoarima)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(1,1,1)(1,1,2)[4]
## Q* = 10.605, df = 3, p-value = 0.01407
##
## Model df: 5. Total lags used: 8
# The residuals aren't like white noise.
# try using other models.
checkresiduals(qauselec_arima0.0.1.0.1.1.4)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,0,1)(0,1,1)[4]
## Q* = 84.595, df = 6, p-value = 4.441e-16
##
## Model df: 2. Total lags used: 8
checkresiduals(qauselec_arima0.1.1.0.1.1.4)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,1,1)(0,1,1)[4]
## Q* = 15.102, df = 6, p-value = 0.01948
##
## Model df: 2. Total lags used: 8
checkresiduals(qauselec_arima0.1.1.0.1.2.4)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(0,1,1)(0,1,2)[4]
## Q* = 13.684, df = 5, p-value = 0.01775
##
## Model df: 3. Total lags used: 8
# The residuals don't resemble white noise regardless of the model. Therefore I'm going to use the best model.
# e. Forecast the next 24 months of data using your preferred model.
fc_qauselec_autoarima <- forecast(
qauselec_autoarima, h = 8
)
autoplot(fc_qauselec_autoarima)
# The forecasts are reasonable.
# f. Compare the forecasts obtained using ets().
fc_qauselec_ets <- forecast(
ets(qauselec), h = 8
)
autoplot(fc_qauselec_ets)
# These forecasts also are reasonable.
14. For the same time series you used in the previous exercise, try using a non-seasonal model applied to the seasonally adjusted data obtained from STL. The stlf() function will make the calculations easy (with method=“arima”). Compare the forecasts with those obtained in the previous exercise. Which do you think is the best approach?
fc_qauselec_stlf <- stlf(
qauselec, lambda = BoxCox.lambda(qauselec),
s.window = 5, robust = TRUE, method = "arima",
h = 8
)
autoplot(fc_qauselec_stlf) +
scale_x_continuous(limits = c(2005, 2012)) +
scale_y_continuous(limits = c(50, 70))
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Warning: Removed 196 row(s) containing missing values (geom_path).
autoplot(fc_qauselec_ets) +
scale_x_continuous(limits = c(2005, 2012)) +
scale_y_continuous(limits = c(50, 70))
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Warning: Removed 196 row(s) containing missing values (geom_path).
autoplot(fc_qauselec_autoarima) +
scale_x_continuous(limits = c(2005, 2012)) +
scale_y_continuous(limits = c(50, 70))
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Warning: Removed 196 row(s) containing missing values (geom_path).
# I don't know which forecasts are best. Forecasts from STL + ARIMA(0, 1, 1) with drift model yielded highest forecasts. It looked like it followed latest trend a lot. ETS(M, A, M) model yielded broadest PI. ARIMA(1, 1, 1)(1, 1, 2)[4] model yielded least varianced forecasts.
# For qauselec data, the choice of model didn't affect forecasts much because they already have strong seasonality and almost non-changing trend.
15. For your retail time series (Exercise 5 above):
# a. develop an appropriate seasonal ARIMA model;
fc_retail_autoarima <- forecast(
auto.arima(retail.ts),
h = 36
)
autoplot(fc_retail_autoarima)
# ARIMA(1, 0, 2)(0, 1, 1)[12] with drift model was chosen.
# b. compare the forecasts with those you obtained in earlier chapters;
# In chapter 3, I used seasonal naive method to forecast.
fc_retail_snaive <- snaive(retail.ts, h = 36)
autoplot(fc_retail_snaive)
# In chapter 7, I thought that Holt-Winters' multiplicative method was best among ets models.
fc_retail_ets <- forecast(
ets(retail.ts, lambda = BoxCox.lambda(retail.ts)),
h = 36
)
autoplot(fc_retail_ets)
# a. Produce a time plot of the sheep population of England and Wales from 1867-1939 (data set sheep).
autoplot(sheep)
# There is decreasing trend, but no definite seasonality.
# b. Assume you decide to fit the following model:
# yt = yt-1 + phi1(yt-1 - yt-2) + phi2(yt-2 - yt-3) + phi3(yt-3 - yt-4) + et
# where et is a white noise series. What sort of ARIMA model is this (i.e., what are p, d, and q)?
# (yt - yt-1) - phi1(yt-1 - yt-2) - phi2(yt-2 - yt-3) - phi3(yt-3 - yt-4) = et
# (1 - B)yt - phi1*B(1- B)yt - phi2*B^2(1- B)yt - phi3*B^3(1- B)yt = et
# (1 - phi1*B - phi2*B^2 - phi3*B^3)(1 - B)yt = et
# It is ARIMA(3, 1, 0) model.
# c. By examining the ACF and PACF of the differenced data, explain why this model is appropriate.
ggtsdisplay(diff(sheep))
sheep.1940 = 1797 + 0.42*(1797 - 1791) -0.20*(1791 - 1627) - 0.30*(1627 - 1665)
sheep.1941 = sheep.1940 + 0.42*(sheep.1940 - 1797) -0.20*(1797 - 1791) - 0.30*(1791 - 1627)
sheep.1942 = sheep.1941 + 0.42*(sheep.1941 - sheep.1940) -0.20*(sheep.1940 - 1797) - 0.30*(1797 - 1791)
c(sheep.1940, sheep.1941, sheep.1942)
## [1] 1778.120 1719.790 1697.268
# e. Now fit the model in R and obtain the forecasts using forecast. How are they different from yours? Why?
fc_sheep_arima.3.1.0 <- forecast(
Arima(sheep, order = c(3, 1, 0)),
h = 3
)
fc_sheep_arima.3.1.0$mean
## Time Series:
## Start = 1940
## End = 1942
## Frequency = 1
## [1] 1777.996 1718.869 1695.985
ar1 <- fc_sheep_arima.3.1.0$model$coef[1]
ar2 <- fc_sheep_arima.3.1.0$model$coef[2]
ar3 <- fc_sheep_arima.3.1.0$model$coef[3]
sheep.1940.new = 1797 + ar1*(1797 - 1791) + ar2*(1791 - 1627) + ar3*(1627 - 1665)
sheep.1941.new = sheep.1940.new + ar1*(sheep.1940.new - 1797) + ar2*(1797 - 1791) + ar3*(1791 - 1627)
sheep.1942.new = sheep.1941.new + ar1*(sheep.1941.new - sheep.1940.new) + ar2*(sheep.1940.new - 1797) + ar3*(1797 - 1791)
c(sheep.1940.new, sheep.1941.new, sheep.1942.new)
## ar1 ar1 ar1
## 1777.996 1718.869 1695.985
# above calculation confirms what I said about the differences.
17.
# a. Plot the annual bituminous coal production in the United States from 1920 to 1968 (data set bicoal).
autoplot(bicoal)
# It looked like there isn't any particular trend or seasonality.
# b. You decide to fit the following model to the series:
# yt = c + phi1*yt-1 + phi2*yt-2 + phi3*yt-3 + phi4*yt-4 + et
# where yt is the coal production in year t and et is a white noise series. What sort of ARIMA model is this (i.e., what are p, d, and q)?
# (1 - phi1*B - phi2*B^2 - phi3*B^3 - phi4*B^4)*yt = c + et
# if mu is the mean of yt,
# c = mu*(1 - phi1*B - phi2*B^2 - phi3*B^3 - phi4*B^4)
# This model is ARIMA(4, 0, 0) or AR(4).
# c. Explain why this model was chosen using the ACF and PACF.
ggAcf(bicoal, lag.max = 36)
ggPacf(bicoal, lag.max = 36)
# ACF plot shows sinusoidally decreasing autocorrelation values. PACF plot shows significant spikes at lag 1 and 4, but none beyond lag 4. Therefore AR(4) model is the appropriate choice.
# d. The last five values of the series are given below.
# Year 1964 1965 1966 1967 1968
# Millions of tons 467 512 534 552 545
# The estimated parameters are c = 162.00, phi1 = 0.83, phi2 = -0.34, phi3 = 0.55 and phi4 = -0.38. Without using the forecast function, calculate forecasts for the next three years (1969-1971).
c = 162.00
phi1 = 0.83
phi2 = -0.34
phi3 = 0.55
phi4 = -0.38
bicoal.1969 <- c + phi1*545 + phi2*552 + phi3*534 + phi4*512
bicoal.1970 <- c + phi1*bicoal.1969 + phi2*545 + phi3*552 + phi4*534
bicoal.1971 <- c + phi1*bicoal.1970 + phi2*bicoal.1969 + phi3*545 + phi4*552
c(bicoal.1969, bicoal.1970, bicoal.1971)
## [1] 525.8100 513.8023 499.6705
# e. Now fit the model in R and obtain the forecasts from the same model. How are they different from yours? Why?
fc_bicoal_ar4 <- forecast(ar(bicoal, 4), h = 3)
fc_bicoal_ar4$mean
## Time Series:
## Start = 1969
## End = 1971
## Frequency = 1
## [1] 526.2057 514.0658 500.0111
# The forecasts from ar function were a little bigger than the calculated forecasts. It also happened because of the small differences of coefficients.
phi1 <- fc_bicoal_ar4$model$ar[1]
phi2 <- fc_bicoal_ar4$model$ar[2]
phi3 <- fc_bicoal_ar4$model$ar[3]
phi4 <- fc_bicoal_ar4$model$ar[4]
c <- fc_bicoal_ar4$model$x.mean*(1 - phi1 - phi2 - phi3 - phi4)
bicoal.1969.new <- c + phi1*545 + phi2*552 + phi3*534 + phi4*512
bicoal.1970.new <- c + phi1*bicoal.1969.new + phi2*545 + phi3*552 + phi4*534
bicoal.1971.new <- c + phi1*bicoal.1970.new + phi2*bicoal.1969.new + phi3*545 + phi4*552
c(bicoal.1969.new, bicoal.1970.new, bicoal.1971.new)
## [1] 526.2057 514.0658 500.0111
# These calculations and results confirm the causation of the differences in forecasts.
